Effects of Neighborhood Environmental Circumstances on Thalamocortical Connectivity

Prepare data

Glasser parcel names

#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) 

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = "orig_parcelname")
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5))
SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382") #colors to use for plotting for each of the 5 quintiles

Environment statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/") #gam output dir, from gam_functions/fit_envGams.R

#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F) 

for(i in 1:length(files)){
  filename <- files[i]
  Rname <- gsub('.{4}$', '', filename)
  x <- read.csv(filename, header = TRUE)
  assign(Rname, x) 
}
rm(x)

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific tract lists

#generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") 
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") 

Associations between Neighborhood and Household SES and Thalamocortical Connectivity

Number of significant environmental effects

##Function to calculate the number and percent of thalamocortical connections showing a significant association with an environmental characteristic
sig.effects <- function(df, env.measure){
  enveffects.df <- df
  enveffects.df$significant <- p.adjust(enveffects.df$GAM.cov.pvalue, method = c("fdr")) < 0.05
  sigeffect.totaln <- sum(enveffects.df$significant)
  sigeffect.percent <- round(sigeffect.totaln/length(enveffects.df$significant)*100, 2)
  print(sprintf("%s thalamocortical connections (%s percent) show a significant association between thalamocortical FA and the %s", sigeffect.totaln, sigeffect.percent, env.measure))
}

PNC

Significant neighborhood-level environment effects (neighborhood SES)

sig.effects(df = envSES_maineffects_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")
## [1] "128 thalamocortical connections (57.4 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"

controlling for parental education

sig.effects(df = envSES_maineffects_educontrolled_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")
## [1] "125 thalamocortical connections (56.05 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"

Percent of positive significant neighborhood-level environment effects

(envSES_maineffects_FA_glasser_pnc %>% mutate(significant = p.adjust(GAM.cov.pvalue, method = c("fdr")) < 0.05) %>% filter(significant == TRUE) %>% filter (GAM.cov.tvalue > 0) %>% nrow() / 128) * 100
## [1] 92.96875

Brain map of neighborhood-level environment effects

envSES_maineffects_FA_glasser_pnc <- envSES_maineffects_FA_glasser_pnc %>% 
  mutate(significant = p.adjust(envSES_maineffects_FA_glasser_pnc$GAM.cov.pvalue, method = c("fdr")) < 0.05)

sigeffects.df <- envSES_maineffects_FA_glasser_pnc %>% filter(significant == TRUE)

sigeffects.positivet.df <- sigeffects.df %>% 
  filter(GAM.cov.tvalue> 0) %>% 
  left_join(glasser.regions, ., by = "tract") %>%
  select(label, tract, GAM.cov.tvalue)
sigeffects.positivet.df$tractlist <- sigeffects.positivet.df$tract %in% tractlist.pnc$tract

ggplot() + 
  geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = GAM.cov.tvalue, colour=I("black"),  size=I(.05))) +  
  theme_void()  + 
  scale_fill_gradientn(colours = c(alpha("#323280", 0.15), "#672975"), limits = c(2, 5.25), oob = squish, na.value = "white") + 
   new_scale_fill() + 
  geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = tractlist, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c("gray89", alpha("white", 0)), na.value = "white") +
   theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/envSES_corticalmap.pdf", device = "pdf", dpi = 500,  width = 4.35, height = 4)

Significant household-level environment effects (average parental years of education)

sig.effects(df = householdSES_maineffects_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")
## [1] "16 thalamocortical connections (7.17 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

controlling for envSES

sig.effects(df = householdSES_maineffects_envSEScontrolled_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

HCPD

Significant household-level environment effects (average parental years of education)

sig.effects(df = householdSES_maineffects_linear_FA_glasser_hcpd, env.measure = "household environment (average parental years of education)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"

Significant household-level environment effects (state-adjusted income-to-needs)

sig.effects(df = INRstate_maineffects_FA_glasser_hcpd, env.measure = "household environment (INR)")
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (INR)"

Thalamocortical Developmental Trajectories Stratified by Neighborhood Environment Factor Score

Trajectories in S-A axis quintiles

#Function to plot predicted developmental trajectories for 10th and 90th percentile envSES factor scores for bins of the S-A axis
plot_trajectories_byenv_SAaxis <- function(quintile, break.points, y.limits){
  
  # Compute average developmental trajectories for low and high factors scores within each S-A axis quintile
  df <- devtrajectories_byenvSES_FA_glasser_pnc %>% merge(.,  S.A.axis, by = c("tract"))
  df$age <- round(df$age, 3)
  df$envSES <- as.factor(df$envSES)
  df <- envSES_maineffects_FA_glasser_pnc %>% select(tract, significant) %>% merge(., df, by = "tract")
  df <- df %>% filter(significant == TRUE)
  
  SAaxismooths.by.envSES <- df %>% group_by(envSES, SA.axis.bin, age) %>% 
                         do(est.mean = mean(.$fitted)) %>% 
                         unnest(cols = c(est.mean)) 

  
  # Plot the trajectories for the requested quintile
  SAtrajectories.byenv.plot <- SAaxismooths.by.envSES %>% filter(SA.axis.bin == quintile) %>%
  ggplot(aes(x = age, y = est.mean, group = interaction(envSES, SA.axis.bin), linetype = envSES)) +
  geom_line(color = SA.colorvector[quintile], linewidth = 0.4) +
  theme_classic() +
  labs(x = "\nAge (years)", y = "\n") +
  scale_y_continuous(breaks = break.points, limits = y.limits) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(SAtrajectories.byenv.plot)  
}
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 1) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[1], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile1.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 1, break.points = c(0.41, 0.42, 0.43, 0.44), y.limits = c(0.405, 0.44))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile1.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 2) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[2], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile2.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 2, break.points = c(0.37, 0.38, 0.39, 0.40), y.limits = c(0.365, 0.40))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile2.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 3) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = "#d4bfe3", na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile3.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 3, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.325, 0.36))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile3.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 4) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[4], na.value = "white")
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry tract                 
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <chr>                 
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY thalamus_L_10pp_autot…
## # ℹ 3 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label                     tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack      L_10pp_ROI     282           4

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile4.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label           geometry tract                 
##   <chr> <chr> <chr> <chr>  <chr> <chr>     <MULTIPOLYGON> <chr>                 
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L_10pp          EMPTY thalamus_L_10pp_autot…
## # ℹ 3 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label                     tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack      L_10pp_ROI     282           4
plot_trajectories_byenv_SAaxis(quintile = 4, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.33, 0.36))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile4.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)
left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 5) %>%
  mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
  ggplot() + 
  geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"),  size=I(.05))) +
  theme_void()  + 
  scale_fill_manual(values = SA.colorvector[5], na.value = "white")
## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile5.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)
## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 5, break.points = c(0.34, 0.35, 0.36, 0.37), y.limits = c(0.34, 0.375))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile5.pdf", device = "pdf", dpi = 500, width = 1.24 , height = 1)

Neighborhood Environment Effects Vary across the Sensorimotor-Association Axis

S-A axis correlation

Spearman’s rho

spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)

cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  spin.data$GAM.cov.tvalue and spin.data$SA.axis
## S = 244280, p-value = 0.0005847
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3010667

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.0056

Correlation plot

ggplot() +
  geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
  geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  labs(x="\nS-A axis", y="Environment effect (t)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/EnvEffect_SArank_correlationplot.pdf", device = "pdf", dpi = 500, width = 1.6, height = 1.3)

Anatomical axis correlations

#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180]  <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$z.mni <- hcp_mmp1.0$z.mni*-1 
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)
#Function to compute the correlation between a environment map and an anatomical axis (x,y,z) as well as test the significance of the difference between the correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(axis){
  df <- merge(spin.data, hcp_mmp1.0, by = "label", sort = F)
  completeobs.n <- sum(!is.na(df$GAM.cov.tvalue))
  
  #S-A axis - environment correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df$SA.axis, df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - environment correlation
  anatomical.axis.cor <- cor(df[,axis], df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and env metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df[,axis], y = df$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df$SA.axis, df[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = completeobs.n, alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

Anterior-posterior axis

ggseg(.data = hcp_mmp1.0, atlas = "glasser", mapping = aes(fill = y.mni)) + paletteer::scale_fill_paletteer_c("harrypotter::lunalovegood")

compute_axis_correlation(axis = "y.mni")
##    axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.06846274         0.38685 0.3010667 0.0005890574

Medial-lateral axis

ggseg(.data = hcp_mmp1.0, atlas = "glasser", mapping = aes(fill = x.mni)) + paletteer::scale_fill_paletteer_c("harrypotter::lunalovegood")

compute_axis_correlation(axis = "x.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.1841123         0.10335 0.3010667    0.3319559

Dorsal-ventral axis

ggseg(.data = hcp_mmp1.0, atlas = "glasser", mapping = aes(fill = z.mni)) + paletteer::scale_fill_paletteer_c("harrypotter::lunalovegood")

compute_axis_correlation(axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.2117744         0.07605 0.3010667    0.3909002

FDR-corrected p-values across cocor comparisons

anatomical.ps <- c(0.0005890574, 0.3319559, 0.3909002)
p.adjust(anatomical.ps, method = c("fdr"))
## [1] 0.001767172 0.390900200 0.390900200
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.3010667)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.06846274)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.2117744)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.1841123)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "D-V", "M-L", "A-P"))

ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) + 
  geom_bar(stat='identity') + 
  labs(x="") +
  labs(y="\nEnvironment Effect Correlation") +
  theme_classic()+
  scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/EnvEffect_anatomicalaxes_barplot.pdf", device = "pdf", dpi = 500, width = 1.55, height = 1.3)

Effect magnitude enrichment testing (S-A axis quintiles)

Mean environment effect t-value by S-A axis bin (empirical)

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.bin) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))
envSES_effects_SAaxisbins
## # A tibble: 5 × 2
##   SA.axis.bin mean.tvalue.significant
##         <dbl>                   <dbl>
## 1           1                    2.70
## 2           2                    3.60
## 3           3                    2.79
## 4           4                    3.91
## 5           5                    4.57

Null distribution and enrichment testing by S-A axis quintile

#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted x inputs 
### set up spun (x) and empirical (y) df
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm.empiricaly <- cbind(y, x.spatialperm) %>% as.data.frame()
colnames(x.spatialperm.empiricaly)[1] <- c("empirical.y")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedx <- sapply(x.spatialperm.empiricaly[2:ncol(x.spatialperm.empiricaly)], #for all permuted bins
                                              function(p){
                                                x.spatialperm.empiricaly %>% group_by(as.character(p)) %>% #group by permuted S-A axis bins
                                                do(mean.tvalue.significant = mean(.$empirical.y, na.rm = TRUE)) %>% #calculate mean t-value
                                                unnest(cols = c("mean.tvalue.significant")) %>% t()}) %>% 
                                                t() %>% as.data.frame()
                                                 
envSES_effects_SAaxisbins_permutedx <- envSES_effects_SAaxisbins_permutedx %>% mutate_if(is.character, as.numeric) #make it numeric
envSES_effects_SAaxisbins_permutedx <-  envSES_effects_SAaxisbins_permutedx[ , !c(TRUE,FALSE) ] 
colnames(envSES_effects_SAaxisbins_permutedx) <- c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue") 

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))

histogram.data <- rbind(envSES_effects_SAaxisbins_permutedx, envSES_effects_SAaxisbins_permutedy)
saveRDS(histogram.data, "/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.rds")
histogram.data <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.rds")
envSES_effects_SAaxisbins_permutedx <- histogram.data[1:10000,]
envSES_effects_SAaxisbins_permutedy <- histogram.data[10001:20000,]
nperm = dim(perm.id.full)[2]
#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(quintile, enrichment_type){
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == quintile) %>% select(mean.tvalue.significant)
  
  #Histogram
  SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382")
  
  null.tvalue.histogram <- histogram.data %>% select(sprintf("SA.bin%s.tvalue", quintile)) %>%
    ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[quintile]) +
    theme_classic() +
    geom_vline(xintercept =  SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
    xlab("Environment effect (t-value)") +
    ylab("Number of nulls") +
    scale_x_continuous(breaks = c(2, 3, 4, 5), limits = c(1.25, 5.25)) +
    theme(
     legend.position = "none", 
     axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
     axis.text.y = element_blank(),
     axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
     axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
     axis.line.x = element_line(linewidth = .2), 
     axis.line.y = element_blank(),
     axis.ticks = element_line(linewidth = .2))
    
  #Permutation-based p
    if(enrichment_type == "smaller"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
    }
  
  if(enrichment_type == "larger"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- p.perm.xy
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
  }
  
  permutation.results <- list(null.tvalue.histogram, perm.p)
  return(permutation.results)
}

S-A quintile 1

enveffects_SAbin_enrichment(quintile = 1, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.0929
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile1_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 2

enveffects_SAbin_enrichment(quintile = 2, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.58425
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile2_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 3

enveffects_SAbin_enrichment(quintile = 3, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.0707
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile3_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 4

enveffects_SAbin_enrichment(quintile = 4, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.1892
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile4_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)

S-A quintile 5

enveffects_SAbin_enrichment(quintile = 5, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.0245
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile5_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)